import numpy as np
import time
from scipy.stats import gmean as geomean

def rp(data, num_perm=100):

	def rank(data):
		return np.argsort(np.argsort(data, axis=0), axis=0)+1

	RP = geomean(rank(data), axis=1)
	order_data = rank(data)

	RP_perm = np.empty([data.shape[0], num_perm]);

	for times in range(num_perm):
	    	new_data = np.copy(data)
	    	for i in range(data.shape[1]):
			new_data[:,i] = np.random.permutation(data[:,i])

		RP_perm[:,times] = geomean(rank(new_data), axis=1)



	temp1 = np.concatenate((RP, RP_perm.flatten()), axis=1)
	temp2 = rank(temp1)
	temp2 = temp2[0:data.shape[0]]
	
	order_temp = np.array([int(np.where(np.sort(temp2)==item)[0]) for item in temp2])
		

	count_perm = np.sort(temp2) - np.array(range(0, data.shape[0]))
	count_perm = count_perm[order_temp]
	pval = count_perm/np.float(num_perm * data.shape[0])

	exp_count = count_perm/np.float(num_perm)
	pfp = exp_count/rank(RP)

	return (RP, pfp, pval)





#def geomean(nums):
#    	return (reduce(lambda x, y: x*y, nums))**(1.0/len(nums))





if __name__=='__main__':
	data=np.random.rand(380, 3)

	start = time.time()
	RP, pfp, pval=rp(data, num_perm=1000)
	end = time.time()
	print end-start

	np.savetxt("testRP.csv", data, delimiter=",")
	np.savetxt("RP.csv", RP, delimiter=",")
	np.savetxt("pfp.csv", pfp, delimiter=",")
